---
title: ""
output:
  pdf_document:
      includes:
          in_header: header.tex 
      keep_tex: yes
  html_notebook: default
  html_document: default
---

```{r setup2, include=FALSE}
rm(list=ls())

options("scipen"=15, "digits"=4)
pkg = list("haven","foreign","ggplot2","plyr","ggthemes","stringr","Hmisc","AER","dplyr","lazyeval","stargazer","knitr","directlabels")
invisible(lapply(pkg, require, character.only = TRUE))
source("./look-up-table.r")

  sigfig <- function(vec, n=2){
    ### function to round values to N significant digits
    # input:   vec       vector of numeric
    #          n         integer is the required sigfig
    # output:  outvec    vector of numeric rounded to N sigfig
    
   if (vec<.0009) {
     formatC(signif(vec,digits=n), digits=n,format="fg", flag="#") 
     } else {
     vec
   }
    
  }      # end of function   sigfig
  
  
print.table <- function(table,caption="test",notes="test",label=NULL,size=2,dig=NULL)
{
  comment          <- list()
  comment$pos      <- list()
  comment$pos[[1]] <- c(nrow(table))
  cols <- ncol(table)
  comment$command  <- c(paste("\\hline \\hline \n ",  # we`ll replace all default hlines with this and the ones below
                              paste0(" \\multicolumn{",cols,"}{p{\\linewidth}}{"),notes,"}  \n  ",
                              sep = ""))
  if (is.null(dig)) {
    suppressWarnings(print(xtable::xtable(table,caption=caption),include.rownames=FALSE,comment=FALSE,caption.placement = 'top',
        add.to.row = comment,
        hline.after = c(-1, 0)))  # indicates rows that will contain hlines (the last one was defined up there)
  } else 
  {
   suppressWarnings(print(xtable::xtable(table,caption=caption,digits=dig),include.rownames=FALSE,comment=FALSE,caption.placement = 'top',
        add.to.row = comment,
        hline.after = c(-1, 0)))   
  }
}


print.table.mc <- function(table,caption="test",notes="test",label=NULL,size=2,head="Outcome: & & \\multicolumn{2}{c}{30D Readmission} & \\multicolumn{2}{c}{30D Mortality} & \\multicolumn{2}{c}{365D Mortality}\\\\",inches=7,sz="\\small")
{
  
  comment          <- list()
  comment$pos      <- list(0,nrow(table))
  #comment$pos[[1]] <- c(nrow(table))
  cols <- ncol(table)
  comment$command  <- c(head,paste("\\hline \\hline \n ",  # we`ll replace all default hlines with this and the ones below
                              paste0(" \\multicolumn{",cols,"}{p{",inches,"in}}{"),notes,"}  \n  ",
                              sep = ""))
  suppressWarnings(print(xtable::xtable(table,caption=caption),include.rownames=FALSE,include.colnames=FALSE,comment=FALSE,caption.placement = 'top',size=sz,
        add.to.row = comment,
        hline.after = c(-1, 0,1,1)))  # indicates rows that will contain hlines (the last one was defined up there)
}



is.extrafont.installed <- function(){
  if(is.element("extrafont", installed.packages()[,1])){
    library(extrafont)
    # probably need something here to run font_import()
    return(T)
  }else{
    warning("Library extrafont installed; using system sans/serif libraries as fallback fonts. 
    To enable full font support, run: 
      install.packages('extrafont') 
      font_import()")
    return(F)
  }
}

base_font_family_tufte <- function(){
  if(is.extrafont.installed()){
    library(extrafont)
    tuftefont <- choose_font(c( "Gill Sans", "GillSans", "Verdana", "serif"), quiet = FALSE)  
  }else{
    tuftefont <- "serif"
  }
  return(tuftefont)
}

theme_tufte_revised <- function(base_size = 11, base_family = base_font_family_tufte(), ticks = TRUE) {

  ret <- theme_bw(base_family = base_family, base_size = base_size) + 
    theme(
          axis.line = element_line(color = 'black'),
          axis.title.x = element_text(vjust = -0.3), 
          axis.title.y = element_text(vjust = 0.8),
          legend.background = element_blank(), 
          legend.key = element_blank(), 
          legend.title = element_text(face="plain"),
          panel.background = element_blank(), 
          panel.border = element_blank(),
          panel.grid = element_blank(),
          plot.background = element_blank(),
          strip.background = element_blank()
    )

  if (!ticks) {
    ret <- ret + theme(axis.ticks = element_blank())
  }

  ret
} 


```


```{r figure 1, echo=FALSE, fig.cap="Density Distribution of Composite Process Score, By Year", message=FALSE, warning=FALSE,cache=TRUE}
load("../data/hospital-compare/Process_Measures_Long.Rdata")
Process.Long %>% mutate(year=lubridate::year(date)) %>% filter(year<=2012) %>% 
  mutate(year=paste(year)) -> xx.p
direct.label(qplot(ABS_PROC_CMP,data=xx.p,colour=year,geom="density"),list("top.points",fontfamily="Gill Sans",cex=.75))+
  theme_tufte_revised()+scale_colour_grey( start = 0, end = 0, na.value = "red")+
  theme(text=element_text(family="Gill Sans"))+
  ylab("Density")+xlab("Composite Process Score")

ranks <- Process.Long %>% mutate(year=lubridate::year(date)) %>% filter(year==2008 | year==2012) %>% select(Provider.ID,ABS_PROC_CMP,year) %>% tbl_df() %>% group_by(Provider.ID,year) %>% filter(row_number()==1) %>%  reshape2::melt(id.vars=c("year","Provider.ID")) %>% reshape2::dcast(Provider.ID~variable+year) %>% na.omit() %>% filter(ABS_PROC_CMP_2012>0 & ABS_PROC_CMP_2008>0) %>% ungroup() %>% arrange(ABS_PROC_CMP_2008) %>%  mutate(rank05 = row_number()) %>%  arrange(ABS_PROC_CMP_2012) %>% mutate(rank12=row_number()) 

#ranks %>% ggplot(aes(x=rank05,y=rank12)) +geom_point(alpha=0.2) +geom_smooth(se=F)


```


```{r load,echo=FALSE,message=FALSE,warning=FALSE,cache=FALSE}
# Summary Functions for Results
# est = function(x) sprintf("%.3f",round(x[2,1],5))
# se = function(x) paste0("(",sprintf("%.3f",round(x[2,2],3)),")")
# estse = function(x) c(est(x),se(x))
# 
# esthr = function(x) sprintf("%.3f",round(x[2:5,1],5))
# sehr = function(x) paste0("(",sprintf("%.3f",round(x[2:5,2],3)),")")
# estsehr = function(x) c(esthr(x),sehr(x))

est = function(x) {
  if (x[2,1] < 0.000999) {
      sigfig(x[2,1])
    } else {
      sprintf("%.3f",round(x[2,1],5))
    }
}

se = function(x)  {
  if (x[2,2]<.00099999) {
    val <- sigfig(x[2,2],1)
  } else {
    val <- sprintf("%.3f",round(x[2,2],3))
  }
  paste0("(",val,")")
  }


estse = function(x) c(est(x),se(x))

esthr = function(x) sprintf("%.3f",round(x[2:5,1],5))
sehr = function(x) paste0("(",sprintf("%.3f",round(x[2:5,2],3)),")")
estsehr = function(x) c(esthr(x),sehr(x))


getn  = function(x,i=4) prettyNum(x[[i]][["n"]],big.mark=",")
stats = function(x,i=2) c(prettyNum(sprintf("%.2f",x[[i]][["stats"]]),big.mark=",",digits=3),prettyNum(x[[i]][["n"]],big.mark=","))

## Measures to Include

QQc = c("Hmort30","Hreadm30","satis","process","Hmort30_abs","Hreadm30_abs","satis_abs","process_abs","Hmort30C","Hreadm30C","satisC","processC")
#QQc = c("Hmort30","Hreadm30","satis","process")
QQprocess.indiv = c("ami2","hf1","hf2","hf3","pn6","inf1","inf3","ami4","pn4","hf4","pn2","pn7","ami1")
QQ = c(QQc,QQprocess.indiv)

####
##
# Non-Discretionary Condition Cohort
##
#
version = "5-0"
cohort = "ndscv220"
setwd("/Volumes/gruber2.nber.org/hospital-quality/analysis")


  
##
# Main Analyses
##
  for (qq in c(QQ)) {
    load(paste0("/Volumes/gruber2.nber.org/hospital-quality/results/",cohort,"-ver-",version,"/","fit-",qq,".Rdata"))
     temp = get(paste0("f", qq))
     
     factor <- 1
     if (grepl("_abs",qq)) factor <- 100
    for (zz in c("death30","readmit30","death365"))
    {
      if (zz=="death30") 
      {
        for(cc in 1:2) 
          {
            temp[[cc]][["iv"]][[zz]][,2] = temp[[cc]][["ses"]][[zz]][,2] 
             temp[[cc]][["ols"]][[zz]][,2] = temp[[cc]][["ses.ols"]][[zz]][,2] 
          }
      } else 
      {
        temp[[2]][["iv"]][[zz]][,2] = temp[[2]][["ses"]][[zz]][,2] 
         temp[[2]][["ols"]][[zz]][,2] = temp[[2]][["ses.ols"]][[zz]][,2] 
      }
    }
 
    temp$fs = cbind(temp[[1]][["fs"]][["readmit30"]] %>% estse(),
         temp[[2]][["fs"]][["readmit30"]] %>% estse()
         )
    
    
    temp$ols = cbind(c(stats(temp)[1],paste0("[",stats(temp)[2],"]")),
                     temp[[1]][["ols"]][["readmit30"]] %>% estse(),
             temp[[2]][["ols"]][["readmit30"]] %>% estse(),
             temp[[1]][["ols"]][["death30"]] %>% estse(),
             temp[[2]][["ols"]][["death30"]] %>% estse(),
             temp[[1]][["ols"]][["death365"]] %>% estse(),
             temp[[2]][["ols"]][["death365"]] %>% estse())
    
    temp$iv = cbind(c(stats(temp)[1],paste0("[",stats(temp)[2],"]")),
                    temp[[1]][["iv"]][["readmit30"]] %>% estse(),
             temp[[2]][["iv"]][["readmit30"]] %>% estse(),
             temp[[1]][["iv"]][["death30"]] %>% estse(),
             temp[[2]][["iv"]][["death30"]] %>% estse(),
             temp[[1]][["iv"]][["death365"]] %>% estse(),
             temp[[2]][["iv"]][["death365"]] %>% estse())
    
    assign(paste0("ndsc.f", qq), temp)
         
  }


##
# RESTAT Sensitivity: No CMS Conditions
##
  for (qq in c(QQ)) {
    load(paste0("/Volumes/gruber2.nber.org/hospital-quality/results/",cohort,"no-ami-pn-hf-ver-",version,"/","fit-",qq,".Rdata"))
     temp = get(paste0("f", qq))
     
     factor <- 1
     if (grepl("_abs",qq)) factor <- 100
    for (zz in c("death30","readmit30","death365"))
    {
      if (zz=="death30") 
      {
        for(cc in 1:2) 
          {
            temp[[cc]][["iv"]][[zz]][,2] = temp[[cc]][["ses"]][[zz]][,2] 
             temp[[cc]][["ols"]][[zz]][,2] = temp[[cc]][["ses.ols"]][[zz]][,2] 
          }
      } else 
      {
        temp[[2]][["iv"]][[zz]][,2] = temp[[2]][["ses"]][[zz]][,2] 
         temp[[2]][["ols"]][[zz]][,2] = temp[[2]][["ses.ols"]][[zz]][,2] 
      }
    }
 
    temp$fs = cbind(temp[[1]][["fs"]][["readmit30"]] %>% estse(),
         temp[[2]][["fs"]][["readmit30"]] %>% estse()
         )
    
    
    temp$ols = cbind(c(stats(temp)[1],paste0("[",stats(temp)[2],"]")),
                     temp[[1]][["ols"]][["readmit30"]] %>% estse(),
             temp[[2]][["ols"]][["readmit30"]] %>% estse(),
             temp[[1]][["ols"]][["death30"]] %>% estse(),
             temp[[2]][["ols"]][["death30"]] %>% estse(),
             temp[[1]][["ols"]][["death365"]] %>% estse(),
             temp[[2]][["ols"]][["death365"]] %>% estse())
    
    temp$iv = cbind(c(stats(temp)[1],paste0("[",stats(temp)[2],"]")),
                    temp[[1]][["iv"]][["readmit30"]] %>% estse(),
             temp[[2]][["iv"]][["readmit30"]] %>% estse(),
             temp[[1]][["iv"]][["death30"]] %>% estse(),
             temp[[2]][["iv"]][["death30"]] %>% estse(),
             temp[[1]][["iv"]][["death365"]] %>% estse(),
             temp[[2]][["iv"]][["death365"]] %>% estse())
    
    assign(paste0("nocms.f", qq), temp)
         
  }


##
# Horserace
##
  load(paste0("/Volumes/gruber2.nber.org/hospital-quality/results/",cohort,"-ver-",version,"/","fit-hr",".Rdata"))
  temp = get(paste0("f", "hr"))
  for (zz in c("death30","readmit30","death365"))
    {
    if (zz=="death30") 
    {
      for(cc in 1:2) 
      {
        temp[[cc]][["iv"]][[zz]][,2] = temp[[cc]][["ses"]][[zz]][,2]
        #temp[[cc]][["ols"]][[zz]][,2] = temp[[cc]][["ses.ols"]][[zz]][,2]
      }
    } else 
    {
       temp[[2]][["iv"]][[zz]][,2] = temp[[2]][["ses"]][[zz]][,2]
       #temp[[2]][["ols"]][[zz]][,2] = temp[[2]][["ses.ols"]][[zz]][,2]
    }
  }
  
  temp$ols = cbind(c(rep("",8)),
                   temp[[1]][["ols"]][["readmit30"]] %>% estsehr(),
             temp[[2]][["ols"]][["readmit30"]] %>% estsehr(),
             temp[[1]][["ols"]][["death30"]] %>% estsehr(),
             temp[[2]][["ols"]][["death30"]] %>% estsehr(),
             temp[[1]][["ols"]][["death365"]] %>% estsehr(),
             temp[[2]][["ols"]][["death365"]] %>% estsehr())
    
    temp$iv = cbind(c(rep("",8)),
                    temp[[1]][["iv"]][["readmit30"]] %>% estsehr(),
             temp[[2]][["iv"]][["readmit30"]] %>% estsehr(),
             temp[[1]][["iv"]][["death30"]] %>% estsehr(),
             temp[[2]][["iv"]][["death30"]] %>% estsehr(),
             temp[[1]][["iv"]][["death365"]] %>% estsehr(),
             temp[[2]][["iv"]][["death365"]] %>% estsehr())
    temp$ols = temp$ols[c(1,5,2,6,3,7,4,8),]
    temp$iv = temp$iv[c(1,5,2,6,3,7,4,8),]
  assign(paste0("ndsc.f", "hr"), temp)


####
##
# CMS Condition Cohort
##
####
version = "4-2"
cohort = "cmsv220"
setwd("/Volumes/gruber2.nber.org/hospital-quality/analysis")

QQcms <- c("Hmort30","Hreadm30","satis","process")

for (qq in QQcms)
{
  load(paste0("/Volumes/gruber2.nber.org/hospital-quality/results/",cohort,"-ver-",version,"/","fit-",qq,".Rdata"))
  temp = get(paste0("f",qq))
  for (zz in c("death30","readmit30","death365"))
  {
    if (zz=="death30") 
    {
      for(cc in 1:2) 
        { 
        temp[[cc]][["iv"]][[zz]][,2] = temp[[cc]][["ses"]][[zz]][,2]
        temp[[cc]][["ols"]][[zz]][,2] = temp[[cc]][["ses.ols"]][[zz]][,2]
        }
    } else 
    {
      temp[[2]][["iv"]][[zz]][,2] = temp[[2]][["ses"]][[zz]][,2]
      temp[[2]][["ols"]][[zz]][,2] = temp[[2]][["ses.ols"]][[zz]][,2]
    }
  }
  
  
      temp$fs = cbind(temp[[1]][["fs"]][["readmit30"]] %>% estse(),
         temp[[2]][["fs"]][["readmit30"]] %>% estse()
         )
    
    temp$ols = cbind(c(stats(temp)[1],paste0("[",stats(temp)[2],"]")),
                     temp[[1]][["ols"]][["readmit30"]] %>% estse(),
             temp[[2]][["ols"]][["readmit30"]] %>% estse(),
             temp[[1]][["ols"]][["death30"]] %>% estse(),
             temp[[2]][["ols"]][["death30"]] %>% estse(),
             temp[[1]][["ols"]][["death365"]] %>% estse(),
             temp[[2]][["ols"]][["death365"]] %>% estse())
    
    temp$iv = cbind(c(stats(temp)[1],paste0("[",stats(temp)[2],"]")),
                    temp[[1]][["iv"]][["readmit30"]] %>% estse(),
             temp[[2]][["iv"]][["readmit30"]] %>% estse(),
             temp[[1]][["iv"]][["death30"]] %>% estse(),
             temp[[2]][["iv"]][["death30"]] %>% estse(),
             temp[[1]][["iv"]][["death365"]] %>% estse(),
             temp[[2]][["iv"]][["death365"]] %>% estse())
    
  assign(paste0("cms.f",qq),temp)
}
# Horserace
load(paste0("/Volumes/gruber2.nber.org/hospital-quality/results/",cohort,"-ver-",version,"/","fit-hr",".Rdata"))
temp = get(paste0("f", "hr"))
for (zz in c("death30","readmit30","death365"))
  {
  if (zz=="death30") 
  {
    for(cc in 1:2) 
      {
      temp[[cc]][["iv"]][[zz]][,2] = temp[[cc]][["ses"]][[zz]][,2]
      temp[[cc]][["ols"]][[zz]][,2] = temp[[cc]][["ses.ols"]][[zz]][,2]
    }
  } else 
  {
     temp[[2]][["iv"]][[zz]][,2] = temp[[2]][["ses"]][[zz]][,2]
     temp[[2]][["ols"]][[zz]][,2] = temp[[2]][["ses.ols"]][[zz]][,2]
  }
}

temp$ols = cbind(c(rep("",8)),temp[[1]][["ols"]][["readmit30"]] %>% estsehr(),
             temp[[2]][["ols"]][["readmit30"]] %>% estsehr(),
             temp[[1]][["ols"]][["death30"]] %>% estsehr(),
             temp[[2]][["ols"]][["death30"]] %>% estsehr(),
             temp[[1]][["ols"]][["death365"]] %>% estsehr(),
             temp[[2]][["ols"]][["death365"]] %>% estsehr())
    
    temp$iv = cbind(c(rep("",8)),temp[[1]][["iv"]][["readmit30"]] %>% estsehr(),
             temp[[2]][["iv"]][["readmit30"]] %>% estsehr(),
             temp[[1]][["iv"]][["death30"]] %>% estsehr(),
             temp[[2]][["iv"]][["death30"]] %>% estsehr(),
             temp[[1]][["iv"]][["death365"]] %>% estsehr(),
             temp[[2]][["iv"]][["death365"]] %>% estsehr())
    temp$ols = temp$ols[c(1,5,2,6,3,7,4,8),]
    temp$iv = temp$iv[c(1,5,2,6,3,7,4,8),]
assign(paste0("cms.f", "hr"), temp)
```


```{r balance,echo=FALSE,results='asis',cache=FALSE,warning=FALSE,message=FALSE}

sigfig <- function(vec, n=3){ 
  ### function to round values to N significant digits
  # input:   vec       vector of numeric
  #          n         integer is the required sigfig  
  # output:  outvec    vector of numeric rounded to N sigfig
  
  formatC(signif(vec,digits=n), digits=n,format="fg", flag="#") 
  
}      # end of function   sigfig


star_from_pvalue <- function(x,sig = c(0.001,0.01,0.05)) {

  n_sig <- length(sig)
  stars <- 1:length(sig) %>%
    purrr::map_chr(~paste0(rep("*",.x),collapse=""))  %>%
    rev()

  out <- stars[findInterval(x, c(0,sig,1))]
  out <- ifelse(is.na(out),"",out)
  return(out)
}
ZZ = "satis"
load(paste0("/Volumes/gruber2.nber.org/hospital-quality/results/ndscv220-ver-5-0/balance-satis.Rdata"))
baltemp <- get(paste0("b",ZZ))
library(purrr)
getcoef <- function(x) x[['rf']][[1]][,1][2] 
getpvalue <- function(x) x[['rf']][[1]][,4][2]
stars <- baltemp %>% map_df(~getpvalue(.x)) %>% t() %>% apply(2,star_from_pvalue) %>% as.vector()
satis.coefs0 <- baltemp %>% map_df(~getcoef(.x)) %>% t() 
satis.coefs <- paste0(sigfig(satis.coefs0,2),stars)

ZZ = "process"
load(paste0("/Volumes/gruber2.nber.org/hospital-quality/results/ndscv220-ver-5-0/balance-",ZZ,".Rdata"))
baltemp <- get(paste0("b",ZZ))
stars <- baltemp %>% map_df(~getpvalue(.x)) %>% t() %>% apply(2,star_from_pvalue) %>% as.vector()
process.coefs0 <- baltemp %>% map_df(~getcoef(.x)) %>% t() 
process.coefs <- paste0(sigfig(process.coefs0,2),stars)

ZZ = "Hmort30"
load(paste0("/Volumes/gruber2.nber.org/hospital-quality/results/ndscv220-ver-5-0/balance-",ZZ,".Rdata"))
baltemp <- get(paste0("b",ZZ))
stars <- baltemp %>% map_df(~getpvalue(.x)) %>% t() %>% apply(2,star_from_pvalue) %>% as.vector()
Hmort30.coefs0 <- baltemp %>% map_df(~getcoef(.x)) %>% t() 
Hmort30.coefs <- paste0(sigfig(Hmort30.coefs0,2),stars)

ZZ = "Hreadm30"
load(paste0("/Volumes/gruber2.nber.org/hospital-quality/results/ndscv220-ver-5-0/balance-",ZZ,".Rdata"))
baltemp <- get(paste0("b",ZZ))
stars <- baltemp %>% map_df(~getpvalue(.x)) %>% t() %>% apply(2,star_from_pvalue) %>% as.vector()
Hreadm30.coefs0 <- baltemp %>% map_df(~getcoef(.x)) %>% t() 
Hreadm30.coefs <- paste0(sigfig(Hreadm30.coefs0,2),stars)

rownames(process.coefs0) = QQ.LUT[rownames(process.coefs0)]
Table1 <- cbind.data.frame(rownames(process.coefs0),rownames(satis.coefs0),process.coefs,satis.coefs,Hmort30.coefs,Hreadm30.coefs)
colnames(Table1) <- c("Measure","Variable","Timely and Effective Care","Patient Satisfaction","30-Day Mortality Rate","30-Day Readmission Rate")
Table1 <- Table1 %>% filter(!is.na(Measure))
Table1.Raw <- Table1
Table1 <- Table1.Raw[-grep("dx",Table1.Raw$Variable), ] %>% select(-Variable)

# version = "4-1"
# setwd("/Volumes/gruber2.nber.org/hospital-quality/analysis")
# 
# load(paste0("/Volumes/gruber2.nber.org/hospital-quality/results/","ndscv220","-ver-",version,"/","balance.Rdata"))
# test = balance[-grep("^dx",rownames(balance)),-5]
# rownames(test) = QQ.LUT[rownames(test)]
# for (c in 1:dim(test)[2]) test[,c]= sprintf("%.3f",as.numeric(test[,c]))
# colnames(test) = paste0("c",1:4)
# test[which(rownames(test)=="Ambulance: Payment"),] = as.character(round(as.numeric(test[which(rownames(test)=="Ambulance: Payment"),]),0))
# 
# #kable(test,cap="Balance-CMS Condition Sample",col.names=c("Q1","Q2","Q3","Q4"))
# Table1 <- cbind.data.frame(rownames(test),test)
# colnames(Table1) <- c("Measure","Quartile 1","Quartile 2","Quartile 3","Quartile 4")
# Table1 <- Table1[-which(rownames(Table1)=="distance"),]
# 

# # Diagnosis Table
# test = as.matrix(balance[grep("^dx",rownames(balance)),-5])
# rownames(test) = QQ.LUT[rownames(test)]
# for (c in 1:dim(test)[2]) test[,c]= sprintf("%.3f",as.numeric(test[,c]))
# colnames(test) = paste0("c",1:4)
# #kable(test,cap="Balance-CMS Condition Sample",col.names=c("Q1","Q2","Q3","Q4"))
# Dx.Table <- cbind.data.frame(rownames(test),test)
# colnames(Dx.Table) <- c("Measure","Quartile 1","Quartile 2","Quartile 3","Quartile 4")

Dx.Table <- Table1.Raw[grep("dx",Table1.Raw$Variable), ] %>% select(-Variable)

```


```{r ,include=FALSE}

# OLS Table: 3 Columns, 1 Panel

blank <- t(as.matrix(rep("",8),nrows=1))

A <- rbind(
  blank,
  cbind(c(QQ.LUT["fprocess"],""),ndsc.fprocess$ols),
  blank,
  cbind(c(QQ.LUT["fsatis"],""),ndsc.fsatis$ols),
  blank,
  cbind(c(QQ.LUT["fHmort30"],""),ndsc.fHmort30$ols),
  blank,
  cbind(c(QQ.LUT["fHreadm30"],""),ndsc.fHreadm30$ols),
  blank
  )

OLS.Table <- rbind(blank,A,blank )[,c(1,2,3,5,7)]
OLS.Table[1,]  = c("Quality Measure","Mean [SD]","(1)","(2)","(3)")


# First Stage
blank <- t(as.matrix(rep("",3),nrows=1))

A <- rbind(
  blank,
  cbind(c(paste0("Ambulance Avg: ",QQ.LUT["fprocess"]),""),ndsc.fprocess$fs),
  blank,
  cbind(c(paste0("Ambulance Avg: ",QQ.LUT["fsatis"]),""),ndsc.fsatis$fs),
  blank,
  cbind(c(paste0("Ambulance Avg: ",QQ.LUT["fHmort30"]),""),ndsc.fHmort30$fs),
  blank,
  cbind(c(paste0("Ambulance Avg: ",QQ.LUT["fHreadm30"]),""),ndsc.fHreadm30$fs),
  blank
  )

FS.Table <- rbind(blank,A,blank )
FS.Table[1,]  = c("Quality Measure Instrument","(1)","(2)")


# 2SLS

blank <- t(as.matrix(rep("",8),nrows=1))

A <- rbind(
  blank,
  cbind(c(QQ.LUT["fprocess"],""),ndsc.fprocess$iv),
  blank,
  cbind(c(QQ.LUT["fsatis"],""),ndsc.fsatis$iv),
  blank,
  cbind(c(QQ.LUT["fHmort30"],""),ndsc.fHmort30$iv),
  blank,
  cbind(c(QQ.LUT["fHreadm30"],""),ndsc.fHreadm30$iv),
  blank
)

IV.Table <- rbind(blank,A,blank )[,c(1,2,3,5,7)]
IV.Table[1,]  = c("Quality Measure","Mean [SD]","(1)","(2)","(3)")


# 2SLS: NO CMS CONDITIONS

blank <- t(as.matrix(rep("",8),nrows=1))

A <- rbind(
  blank,
  cbind(c(QQ.LUT["fprocess"],""),nocms.fprocess$iv),
  blank,
  cbind(c(QQ.LUT["fsatis"],""),nocms.fsatis$iv),
  blank,
  cbind(c(QQ.LUT["fHmort30"],""),nocms.fHmort30$iv),
  blank,
  cbind(c(QQ.LUT["fHreadm30"],""),nocms.fHreadm30$iv),
  blank
)

IV.Table.NoCMS <- rbind(blank,A,blank )[,c(1,2,3,5,7)]
IV.Table.NoCMS[1,]  = c("Quality Measure","Mean [SD]","(1)","(2)","(3)")

blank <- t(as.matrix(rep("",8),nrows=1))

A <- rbind(
  blank,
  cbind(c(QQ.LUT["fprocess"],""),ndsc.fprocess$iv),
  blank,
  cbind(c(QQ.LUT["fsatis"],""),ndsc.fsatis$iv),
  blank,
  cbind(c(QQ.LUT["fHmort30"],""),ndsc.fHmort30$iv),
  blank,
  cbind(c(QQ.LUT["fHreadm30"],""),ndsc.fHreadm30$iv),
  blank
)
A <- A[,c(1,2,4,6,8)]
B <- rbind(
  blank,
  cbind(c(QQ.LUT["fprocess"],"",QQ.LUT["fsatis"],"",QQ.LUT["fHmort30"],"",QQ.LUT["fHreadm30"],""),ndsc.fhr$iv),
  blank)[,c(1,2,4,6,8)]

C <- rbind(
  blank,
  cbind(c(QQ.LUT["fprocess"],""),cms.fprocess$iv),
  blank,
  cbind(c(QQ.LUT["fsatis"],""),cms.fsatis$iv),
  blank,
  cbind(c(QQ.LUT["fHmort30"],""),cms.fHmort30$iv),
  blank,
  cbind(c(QQ.LUT["fHreadm30"],""),cms.fHreadm30$iv),
  blank
)
C <- C[,c(1,2,4,6,8)]
blank <- t(as.matrix(C[1,]))
C[1,1] = "Panel C. CMS Condtion Sample (AMI,PN,HF)"
A[1,1] = "Panel A. Add Comorbidity Controls"
B[1,1] = "Panel B. Horse-Race"
B[,2] = A[c(1,2,3,5,6,8,9,11,12,13),2]
IV.Table.2 <- rbind(blank,A,blank,B,blank,C)

# 2SLS : Unstandardized Measures

blank <- t(as.matrix(rep("",8),nrows=1))

A <- rbind(
  blank,
  cbind(c(QQ.LUT["fprocess"],""),ndsc.fprocess_abs$iv),
  blank,
  cbind(c(QQ.LUT["fsatis"],""),ndsc.fsatis_abs$iv),
  blank,
  cbind(c(QQ.LUT["fHmort30"],""),ndsc.fHmort30_abs$iv),
  blank,
  cbind(c(QQ.LUT["fHreadm30"],""),ndsc.fHreadm30_abs$iv),
  blank
)

IV.Table.Abs <- rbind(blank,A,blank )[,c(1,2,3,5,7)]
IV.Table.Abs[1,]  = c("Quality Measure","Mean [SD]","(1)","(2)","(3)")


# 2SLS : Concurrent Quality

blank <- t(as.matrix(rep("",8),nrows=1))

A <- rbind(
  blank,
  cbind(c(QQ.LUT["fprocessC"],""),ndsc.fprocessC$iv),
  blank,
  cbind(c(QQ.LUT["fsatisC"],""),ndsc.fsatisC$iv),
  blank,
  cbind(c(QQ.LUT["fHmort30C"],""),ndsc.fHmort30C$iv),
  blank,
  cbind(c(QQ.LUT["fHreadm30C"],""),ndsc.fHreadm30C$iv),
  blank
)

IV.Table.Concurrent <- rbind(blank,A,blank )[,c(1,2,3,5,7)]
IV.Table.Concurrent[1,]  = c("Quality Measure","Mean [SD]","(1)","(2)","(3)")


# 2SLS : Process measures

blank <- t(as.matrix(rep("",8),nrows=1))

A <- rbind(
  blank,
  cbind(c(QQ.LUT["fami4"],""),ndsc.fami4$iv),
  blank,
  cbind(c(QQ.LUT["fpn4"],""),ndsc.fpn4$iv),
  blank,
  cbind(c(QQ.LUT["fhf4"],""),ndsc.fhf4$iv),
  blank
)

Smoke.Table <- rbind(blank,A,blank )[,c(1,2,3,5,7)]
Smoke.Table[1,]  = c("Quality Measure","Mean [SD]","(1)","(2)","(3)")
```
\blandscape

```{r,echo=FALSE,results='asis'}
print.table(Table1,dig=3,caption="Balance of Patient Characteristics by Quality Instrument",notes=c("Table shows balance across controls used in regressions. Balance is assessed using pairwise regressions of each characteristic on each instrument based on the specification in Equation 5. The reported estimates (rows) show the coefficient on the instrument from each of these regressions fit separately for each instrument (columns).   Sample size = 546,700.  ***(**)(*) = significant at 0.001 (0.01) (0.05) level."))
```




```{r,echo=FALSE,results='asis'}
print.table.mc(table=OLS.Table,size=1,caption= "OLS Results",notes="Each cell reports ordinary least squares (OLS) coefficient estimates for a separate regression.  Quality Measures have been demeaned and standardized by 2 standard deviations so they can be interpreted like binary (low-to-high) indicators. The underlying mean and standard deviation of each quality measure is provided in the first column to facilitate interpretation on the original scale. Outcome means: 30D Readmission = 15.0\\%, 30D Mortality = 17.0\\%, 365D Mortality = 37.2\\%. Sample sizes: 546,700 (30-Day outcomes), 451,503 (1-Year Mortality). All models include patient demographic and ambulance controls as listed in Table 1, as well as the diagnosis controls as listed in Table A2. Models also include ZIP code-patient origin fixed effects. Standard errors, clustered at Health Service Area (HSA) level, are reported in parentheses.",
head="Outcome: & & \\multicolumn{1}{c}{30D Readmission} & \\multicolumn{1}{c}{30D Mortality} & \\multicolumn{1}{c}{365D Mortality}\\\\")

print.table.mc(table=FS.Table,size=1,caption= "First Stage Results",notes="Each cell reflects a separate first-stage regression of the ambulance instrument on the quality measure. Quality Measures have been demeaned and standardized by 2 standard deviations so they can be interpreted like binary (low-to-high) indicators. The underlying mean and standard deviation of each quality measure is provided in the first column to facilitate interpretation on the original scale. Outcome means: 30D Readmission = 15.0\\%, 30D Mortality = 17.0\\%, 365D Mortality = 37.2\\%. Sample sizes: 546,700 (30-Day outcomes), 451,503 (1-Year Mortality). All models include patient demographic and ambulance controls as listed in Table 1, as well as the diagnosis controls as listed in Table A2. Models also include ZIP code-patient origin fixed effects. Comorbidity controls are listed in Table 1. Standard errors, clustered at Health Service Area (HSA) level, are reported in parentheses.",
head="&  \\multicolumn{1}{c}{No Comorbidity Controls} & \\multicolumn{1}{c}{With Comorbidity Controls} \\\\",inches=7)

print.table.mc(table=IV.Table,size=1,caption= "2SLS Results",notes="Each cell reports two-stage least squares (2SLS) coefficient estimates for a separate regression.  Quality Measures have been demeaned and standardized by 2 standard deviations so they can be interpreted like binary (low-to-high) indicators. The underlying mean and standard deviation of each quality measure is provided the first column to facilitate interpretation on the original scale. Outcome means: 30D Readmission = 15.0\\%, 30D Mortality = 17.0\\%, 365D Mortality = 37.2\\%. Sample sizes: 546,700 (30-Day outcomes), 451,503 (1-Year Mortality). All models include patient demographic and ambulance controls as listed in Table 1, as well as the diagnosis controls as listed in Table A2. Models also include ZIP code-patient origin fixed effects. Standard errors, clustered at Health Service Area (HSA) level, are reported in parentheses.",
               head="Outcome: & & \\multicolumn{1}{c}{30D Readmission} & \\multicolumn{1}{c}{30D Mortality} & \\multicolumn{1}{c}{365D Mortality}\\\\")
```

\elandscape

```{r,results='asis',echo=FALSE}
print.table.mc(table=IV.Table.2,size=1,caption= "2SLS Results",notes="In Panel A, each cell reports two-stage least squares (2SLS) coefficient estimates for a separate regression. In Panel B, each column reports 2SLS estimates for a 'horse race' specificiation that includes all quality measures in a single regression. In Panel C, each cell reports two-stage least squares (2SLS) coefficient estimates for a separate regression using a sample of only Acute Myocardial Infarction (AMI), Pneumonia (PN) and Heart Failure (HF) patients. All models include patient demographic, comorbidity and ambulance controls as listed in Table 1, as well as the diagnosis controls as listed in Table A2. Models also include ZIP code-patient origin fixed effects. Quality Measures have been demeaned and standardized by 2 standard deviations so they can be interpreted like binary (low-to-high) indicators. The underlying mean and standard deviation of each quality measure is provided the first column to facilitate interpretation on the original scale. Outcome means for non-discretionary condition sample: 30D Readmission = 15.0\\%, 30D Mortality = 17.0\\%, 365D Mortality = 37.2\\%. Sample sizes for non-discretionary condition sample: 546,700 (30-Day outcomes), 451,503 (1-Year Mortality).  Outcome means for CMS condition sample: 30D Readmission = 19.0\\%, 30D Mortality = 15.8\\%, 365D Mortality = 41.0\\%. Sample sizes for CMS condition sample: 171,246 (30-Day outcomes), 142,424 (1-Year Mortality). Standard errors, clustered at Health Service Area (HSA) level, are reported in parentheses. ",head="Outcome: & Mean [SD] & \\multicolumn{1}{c}{30D Readmission} & \\multicolumn{1}{c}{30D Mortality} & \\multicolumn{1}{c}{365D Mortality}\\\\")
```

\blandscape


```{r, echo=FALSE, message = FALSE, warning = FALSE, cache=FALSE,results='asis'}
qq <- read_dta("../data/aha-data-112116.dta") %>% filter(year==2012) %>% select(provider.id,
    satis = abs_hcahps_cmp_yr1,
  process = abs_proc_cmp_yr1,
  Hmort30 = abs_mort_cmp_yr1,
  Hreadm30 = abs_readm_cmp_yr1,
  AMImort30 = mort_30_ami_yr1,
  PNmort30 = mort_30_pn_yr1,
  HFmort30 = mort_30_hf_yr1,
  AMIreadm30 = readm_30_ami_yr1,
  PNreadm30 = readm_30_pn_yr1,
  HFreadm30 = readm_30_hf_yr1,
  
  volume = hospbd,
  teach = teach,
  forpr = forpr,
  nonpr = nonpr,
  gov = gov,
  lowpr = low.profit,
  hipr = high.profit,
  coth = coth,
  diag_year = year
  )
load("/Volumes/gruber2.nber.org/hospital-quality/data/hospital-composite-measures-ndscv220-ver-4-1.Rdata")
composite <- AAcomp %>% rename(provider.id = prvnumgrp)
load("/Volumes/gruber2.nber.org/hospital-quality/data/hospital-sample-cmsv220-ver-4-0.Rdata")
qq <- qq %>% mutate(prvnumgrp=provider.id) %>% inner_join(AAhosps,"prvnumgrp") %>% mutate(volume.st = (volume-mean(volume))/sd(volume)) %>% inner_join(composite,c("provider.id","diag_year"))
fit.teach <- with(qq,lm(process~teach))
fit.volume <- with(qq,lm(process~volume.st))
fit.own <- with(qq,lm(process~forpr+gov))
fit.all.process <- with(qq,lm(process~teach+volume.st+forpr+gov))
fit.all.satis <- with(qq,lm(satis~teach+volume.st+forpr+gov))
fit.all.mort <- with(qq,lm(Hmort30~teach+volume.st+forpr+gov))
fit.all.readm <- with(qq,lm(Hreadm30~teach+volume.st+forpr+gov))

# stargazer(list(fit.all.process,fit.all.satis,fit.all.mort,fit.all.readm),title="Association Between  Composite Quality Scores and Structural Hospital Characteristics",covariate.labels=c("Teaching Hospital","Hospital Size (Standardized)","Ownership: Investor-Owned","Ownership: Public"),
#           omit.stat = c("adj.rsq","rsq","f","ser"),
#           dep.var.caption=c(""),
#           dep.var.labels=c(""),
#            dep.var.labels.include=FALSE,
#           model.names=FALSE,
#           keep.stat=c("sigma2"),
#           notes.align="l",
#           star.char=NULL,
#           font.size="small",
#           star.cutoffs= NULL,
#           table.placement= "!h",header=FALSE,notes=c("Standard Errors in parentheses."),report=c("vcs"),
#           column.labels=c("Timely and Effective Care","Patient Satisfaction","30D Mortality","30D Readmission")
#           )

```


<!-- \begin{table}[!h] \centering  -->
<!--   \caption{Association Between  Composite Quality Scores and Structural Hospital Characteristics, 2012 Values}  -->
<!--   \label{}  -->
<!-- \small  -->
<!-- \begin{tabular}{@{\extracolsep{5pt}}lcccc}  -->
<!-- \\[-1.8ex]\hline  -->
<!-- \hline \\[-1.8ex]  -->
<!--  & Timely and Effective Care & Patient Satisfaction & 30D Mortality & 30D Readmission \\  -->
<!-- \\[-1.8ex] & (1) & (2) & (3) & (4)\\  -->
<!-- \hline \\[-1.8ex]  -->
<!--  Teaching Hospital & 0.427 & $-$0.027 & $-$0.140 & 0.014 \\  -->
<!--   & (0.166) & (0.239) & (0.066) & (0.076) \\  -->
<!--   & & & & \\  -->
<!--  Hospital Beds (Standardized) & 0.353 & $-$0.518 & $-$0.214 & 0.175 \\  -->
<!--   & (0.078) & (0.112) & (0.031) & (0.035) \\  -->
<!--   & & & & \\  -->
<!--  Ownership: Investor-Owned & 1.013 & $-$2.105 & 0.263 & 0.356 \\  -->
<!--   & (0.175) & (0.251) & (0.070) & (0.080) \\  -->
<!--   & & & & \\  -->
<!--  Ownership: Public & $-$1.590 & 0.002 & 0.587 & 0.016 \\  -->
<!--   & (0.211) & (0.303) & (0.084) & (0.096) \\  -->
<!--   & & & & \\  -->
<!--  Constant & 96.540 & 69.790 & 12.880 & 20.930 \\  -->
<!--   & (0.099) & (0.143) & (0.040) & (0.045) \\  -->
<!--   & & & & \\  -->
<!-- \hline \\[-1.8ex]  -->
<!-- \hline  -->
<!-- \hline \\[-1.8ex]  -->
<!-- \multicolumn{4}{p{7in}}{Table shows results of a ``horse race'' OLS regression of each quality measure on hospital characteristics. The dependent variable in each regression is expressed in terms of standard deviations. Sample size = 2,116 general acute care hospitals.} \\  -->
<!-- \end{tabular}  -->
<!-- \end{table}  -->



\setcounter{table}{0}
\renewcommand*\thetable{A\arabic{table}}
\renewcommand{\thefigure}{A\arabic{figure}}

```{r,echo=FALSE,results='asis'}

blank <- t(as.matrix(rep("",1),nrows=1))

A <- rbind(
  blank,
  "Heart Failure Patients Given ACE Inhibitor or ARB for Left Ventricular Systolic Dysfunction (LVSD)",
  "AMI Patients Given Aspirin at Discharge",
  "Heart Failure Patients Given Assessment of Left Ventricular Function (LVF)",
  "Heart Failure Patients Given Discharge Instructions",
  "Pneumonia Patients Given the Most Appropriate Initial Antibiotic(s)",
  "Surgery Patients Who Received Preventative Antibiotic(s) One Hour Before Incision",
  "Surgery Patients Whose Preventative Antibiotic(s) are Stopped Within 24 hours After Surgery",
  blank
)

B <- rbind(
  blank,
  "Doctors always communicated well",
  "Nurses always communicated well",
  "Pain was always well controlled",
  "Patients always received help as soon as they wanted – Patients who gave a rating of 9 or 10 (high)",
  "Room was always clean",  
  "Room always quiet at night",
  "Staff always explained",
  "Yes, patients would definitely recommend the hospital ",
  "Yes, staff did give patients this information",
  blank
)


C <- rbind(
  blank,
  "AMI 30-Day Mortality Rate",
  "Pneumonia 30-Day Mortality Rate",
  "Heart Failure 30-Day Mortality Rate",blank
)

D <- rbind(
  blank,
  "AMI 30-Day Readmission Rate",
  "Pneumonia 30-Day Readmission Rate",
  "Heart Failure 30-Day Readmission Rate",blank
)


blank <- t(as.matrix(C[1,]))
Measure.Table <- data.frame(rbind(blank,A,blank,B,blank,C,blank,D))
names(Measure.Table) = c("Measure")
Measure.Table$Domain = ""
Measure.Table <- as.matrix(Measure.Table[,c(2,1)])
Measure.Table[2,1] = "Timely and Effective Care"
Measure.Table[12,1] = "Patient Experience of Care"
Measure.Table[24,1] = "30-Day Mortality Rates"
Measure.Table[30,1] = "30-Day Readmission Rates"


print.table(Measure.Table,caption="CMS Quality Measures Used to Construct Composite Domain Scores",notes="Source: Archived Hospital Compare website data.  The process and outcome measures are defined for specific patient conditions (e.g., AMI, pneumonia, heart failure) whereas the HCAHPS patient satisfaction surveys are distributed to a wider range of patients. ")
```

```{r,echo=FALSE,results='asis'}
print.table(Dx.Table[,c(1,2,3)],dig=3,caption="Balance: Non-Deferrable Condition Sample",notes=c("Table shows balance across principal diagnosis controls used in regressions. Balance is assessed using pairwise regressions of each characteristic on each instrument based on the specification in Equation 5. The reported estimates (rows) show the coefficient on the instrument from each of these regressions fit separately for each instrument (columns).   Sample size = 546,700"))

```
```{r,echo=FALSE,results='asis'}
print.table(Dx.Table[,c(1,4,5)],dig=3,caption="Balance: Non-Deferrable Condition Sample",notes=c("Table shows balance across principal diagnosis controls used in regressions. Balance is assessed using pairwise regressions of each characteristic on each instrument based on the specification in Equation 5. The reported estimates (rows) show the coefficient on the instrument from each of these regressions fit separately for each instrument (columns).   Sample size = 546,700"))

```


```{r,echo=FALSE,results='asis'}
print.table.mc(table=IV.Table.Abs,size=1,caption= "2SLS Results - Unstandardized Measures",notes="Each cell reports two-stage least squares (2SLS) coefficient estimates for a separate regression.  Outcome means: 30D Readmission = 15.0\\%, 30D Mortality = 17.0\\%, 365D Mortality = 37.2\\%. Sample sizes: 546,700 (30-Day outcomes), 451,503 (1-Year Mortality). All models include patient demographic and ambulance controls as listed in Table 1, as well as the diagnosis controls as listed in Table A2. Models also include ZIP code-patient origin fixed effects. Standard errors, clustered at Health Service Area (HSA) level, are reported in parentheses.",
               head="Outcome: & & \\multicolumn{1}{c}{30D Readmission} & \\multicolumn{1}{c}{30D Mortality} & \\multicolumn{1}{c}{365D Mortality}\\\\")
```


<!-- # ```{r, echo=FALSE, message = FALSE, warning = FALSE, cache=FALSE,results='asis'} -->
<!-- # print.table.mc(table=Smoke.Table,size=1,caption= "2SLS Results: Alternative Process Measures",notes="Each cell reports two-stage least squares (2SLS) coefficient estimates for a separate regression.  Quality Measures have been demeaned and standardized by 2 standard deviations so they can be interpreted like binary (low-to-high) indicators. The underlying mean and standard deviation of each quality measure is provided the first column to facilitate interpretation on the original scale. Outcome means: 30D Readmission = 15.0\\%, 30D Mortality = 17.0\\%, 365D Mortality = 37.2\\%. Sample sizes: 546,700 (30-Day outcomes), 451,503 (1-Year Mortality). All models include patient demographic and ambulance controls as listed in Table 1, as well as the diagnosis controls as listed in Table A2. Models also include ZIP code-patient origin fixed effects. Standard errors, clustered at Health Service Area (HSA) level, are reported in parentheses.", -->
<!-- #                head="Outcome: & & \\multicolumn{1}{c}{30D Readmission} & \\multicolumn{1}{c}{30D Mortality} & \\multicolumn{1}{c}{365D Mortality}\\\\") -->
<!-- # ``` -->


```{r,echo=FALSE,results='asis'}
print.table.mc(table=IV.Table.Concurrent,size=1,caption= "2SLS Results - Concurrent Quality Measures",notes="This table replicates the results from Table 4 but using measures of quality with measurement that runs concurrent to the patient data.  For example, the process quality measure used for patients admitted in 2008 reflects the values published in 2009, since those values were measured based on 2008 admissions. Each cell reports two-stage least squares (2SLS) coefficient estimates for a separate regression.  Outcome means: 30D Readmission = 15.0\\%, 30D Mortality = 17.0\\%, 365D Mortality = 37.2\\%. Sample sizes: 546,700 (30-Day outcomes), 451,503 (1-Year Mortality). All models include patient demographic and ambulance controls as listed in Table 1, as well as the diagnosis controls as listed in Table A2. Models also include ZIP code-patient origin fixed effects. Standard errors, clustered at Health Service Area (HSA) level, are reported in parentheses.",
               head="Outcome: & & \\multicolumn{1}{c}{30D Readmission} & \\multicolumn{1}{c}{30D Mortality} & \\multicolumn{1}{c}{365D Mortality}\\\\")
```




```{r,echo=FALSE,results='asis'}
print.table.mc(table=IV.Table.NoCMS,size=1,caption= "2SLS Results - Excluding AMI, and PN Patients",notes="This table replicates the results from Table 4 but using nondeferrable conditions that do not include AMI, HF, and PN patients. Each cell reports two-stage least squares (2SLS) coefficient estimates for a separate regression.  Outcome means: 30D Readmission = 14.5\\%, 30D Mortality = 17.2\\%, 365D Mortality = 37.1\\%. Sample sizes: 438,911 (30-Day outcomes), 361,958 (1-Year Mortality). All models include patient demographic and ambulance controls as listed in Table 1, as well as the diagnosis controls as listed in Table A2. Models also include ZIP code-patient origin fixed effects. Standard errors, clustered at Health Service Area (HSA) level, are reported in parentheses.",
               head="Outcome: & & \\multicolumn{1}{c}{30D Readmission} & \\multicolumn{1}{c}{30D Mortality} & \\multicolumn{1}{c}{365D Mortality}\\\\")
```

\elandscape


